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Abstract —Three-dimensional x-ray CT image reconstruction 
in baggage scanning in security applications is an important 
research field. The variety of materials to be reconstructed is 
broader than medical x-ray imaging. Presence of high atten¬ 
uating materials such as metal may cause artifacts if analyti¬ 
cal reconstruction methods are used. Statistical modeling and 
the resultant iterative algorithms are known to reduce these 
artifacts and present good quantitative accuracy in estimates 
of linear attenuation coefficients. However, iterative algorithms 
may require computations in order to achieve quantitatively 
accurate results. For the case of baggage scanning, in order 
to provide fast accurate inspection throughput, they must be 
accelerated drastically. There are many approaches proposed 
in the literature to increase speed of convergence. This paper 
presents a new method that estimates the wavelet coefficients of 
the images in the discrete wavelet transform domain instead of 
the image space itself. Initially, surrogate functions are created 
around approximation coefficients only. As the iterations proceed, 
the wavelet tree on which the updates are made is expanded 
based on a criterion and detail coefficients at each level are 
updated and the tree is expanded this way. For example, in 
the smooth regions of the image the detail coefficients are not 
updated while the coefficients that represent the high-frequency 
component around edges are being updated, thus saving time by 
focusing computations where they are needed. This approach is 
implemented on real data from a SureScan™ jclOOO Explosive 
Detection SystenQand compared to straightforward implementa¬ 
tion of the unregularized alternating minimization of O’Sullivan 
and Benac Q. 

I. Summary 

X-ray CT image reconstruction algorithms often are de¬ 
signed to trade-off a data fit term and an image roughness term. 
Typical examples of data fit terms are squared error 0 and 
log-likelihood. Typical examples of image roughness terms 
are total variation, Gauss-Markov random field priors 0, and 
Huber class roughness penalties. A complementary approach 
for penalizing roughness is to represent the image using a 
wavelet (or other multiresolution) expansion, directly estimate 
the wavelet coefficients, and introduce a penalty (often Li) on 
the wavelet coefficients 0 , 0 . Multigrid approach 0 - 11 ) 

^ SureScan™ is a trademark of the SureScan Corporation. 


has been shown to be successful in some image reconstruction 
methods in terms of achieving faster convergence speed. The 
idea is to move through different grid levels over time. At any 
grid level, the voxel size is the same throughout the image 
domain. The approach we describe in this paper is closest to 
this with the following differences: 

• Our data fit term is a Poisson log-likelihood with mean 
determined by Beer’s law 

• We use an alternating minimization framework to derive 
an algorithm that is guaranteed to decrease the cost 
function at every iteration 

• We extend the prior alternating minimization framework 
to point spread functions with negative values 

• We update only a subset of wavelet coefficients, con¬ 
strained to be on a tree, thereby decreasing the computa¬ 
tional complexity per iteration relative to fully updating 
the image 

• We adaptively update the tree defining which wavelet 
coefficients are updated at each iteration 

• Our wavelet tree structure results in image domain rep¬ 
resentation of voxels with different sizes 

• We incorporate an adaptive threshold allowing the com¬ 
putation of a sequence of images of increasing resolution, 
with increasing roughness and increasing log-likelihood 

The result is a fast, adaptive, iterative image reconstruction 
algorithm. 

II. Introduction 

The problem to be minimized for x-ray CT in this paper is 
penalized likelihood estimation with a Poisson log-likelihood 
data fitting term and a regularization term, optimized over 
image G It is shown in Q that maximization 

of the Poisson log-likelihood term is equivalent to mini¬ 
mizing the I-divergenc^ between the transmission data d 
and the estimated mean q(/x) G , where q{/ji){y) = 

^I-divergence between two vectors p,q E is defined as I{p\\q) = 
J2iPilog{^) -Pi + qi. 


Io{y) eyiY>{—^^h{y\x)iJi{x)), /o(.) is the incident photon 
count vector, h{y\x) is an element of the system matrix 
H G that represents the length of the intersection 

between the ray path of index y G and voxel of index 
X G . Then this penalized likelihood estimation problem 
can be formulated as GD 


Mpml =argmin/(d||q(p)) + Ai?(p), (1) 

A ^>0 

where R{iJi) is a regularization term selected as a roughness 
penalty and A > 0 is the parameter that controls the level of 
roughness imposed on the image. Also, it is important to note 
that the non-negativity constraint on fi is due to the nature of 
linear attenuation coefficients of materials. Since there is no 
closed form solution to this problem, we solve it iteratively. 
At each iteration, a surrogate function that approximates 
the original objective function is minimized, which in turn 
decreases the original objective function. In our recent work 
|TT|, we generalized the formulation of surrogate functions in 
I y for data fitting term to the regularization term. The idea is 
to use Jensen’s inequality to decouple the objective function 
and form many one-parameter convex functions, minimize 
them, and iterate. 

Assume that there exists a discrete wavelet inverse transform 
matrix ft G that is non-singular. Then, the image fi 

can be represented as 


M (2) 

where j3 is the vector of wavelet coefficients. The problem in 
this paper can then be written as 

l^*PML = argmm/(d||q(f2/3)) + Ai?(f2/3) (3) 
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subject to r2/3 > 0 

Below, the derivation of the surrogate functions for the data 
fitting term is shown. A similar approach yields surrogate 
functions for the regularization term as well. 

The I-divergence term can be written as 

I(d\\q{^^|3)) = '^d{y)'^h{y\x)'^u{z\y)l3{z) 

y X z 

+ (“ X! Hy\x)'^^ix\z)i3(z)) 

y X z 

-{-constantly). (4) 


follows. 
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+ X] d{y) exp (- 4‘{y\z)W{z) - /3(«)) 
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( 5 ) 


where 

b(z) = '^d{y)(j)(y\z), (6) 

y 

the convex decomposition lemma Q is used for r{z\y) > 0, 

chosen as 


(myMl 
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+ yy const{z') (8) 

z'tz. 


Adding the constant term in I-divergence, we define our 
surrogate function. 


4,(d||q;/3,/3)= y] b{z)P{z) 

+ const{z') + const(y) (9) 

z'^Z^ 


It is clear to see that this is a one-parameter convex function 
over each ^(z) and the gradient with respect to ^[z) is given 
as: 


For simplicity, define the matrix ^ = Hft, where 

(l){y\z) is the system matrix element between ray path 
of index y and wavelet coefficient of index 2 ; G 
Z^. Assume that there exists a known estimate /3 
and q{y) = Io{y)&zp{-Y.xKv\x)Y.z‘^i^\^)Pi^)) = 

/o(y) exp(—(/)(y|2;)/3(2;)). The terms in the I-divergence 
that depend on (3 are used to construct surrogate functions as 


4(d||q;/3,/3) = b{z)-h+{z)eM-Zmz)-P{z))) 

- S_(2)exp(2o(/3(«)-^(«))) (10) 

^Zs represents a subset of the wavelet domain to be chosen for update. 
In our approach, we choose it in a way that every voxel in image domain is 
represented at any iteration, possibly with different numbers of coefficients. 
This subset can be fixed or be varied over iterations. 









where 


hiz) = 

qiy)(j)iy\z), 

(11) 

yMy\z)>o 



b-iz)= E 

4{y)<t>{y\z)- 

(12) 


y,0(y|z)<O 


The first-order necessary condition for a minimizer is to find 
the ^{z) for which the gradient is zero, which has a closed 
form solution. The algorithm is shown below. 



Algorithm 1 Unregularized Wavelet AM Algorithm 

Inputs: Iq, H, f2, zi^horj = 0,1,(J - 1) 

Precompute b(z) = Y,y d(y)(l>{y\z) 
for j = 0,1,(J — 1) do 

g(j)(j/) = exp(- (I>iy\z)l3^^\z)) 

=maxj,E^gzO) \(l>iy\z)\ 

(i) ^ 

for every 2; G Zs do 

b+Hz) = T,y,^(y\z)>0^iy)<l>iy\^) 

b^Hz) = T,y,^(y\z)<0^iy)<l>iy\^) 

^(j+i)(2;) = p(^z) where 

biz) - b^i\z) exp(-Z«(/3(^) - /?0)(^))) 

-b^l\z)eMZo\^iz)-P^^Hz))) = 0 

end for 
end for 


Fig. 1: Objective function values vs. time for AM and Wavelet 
AM. 



Fig. 2: Image reconstructed with unregularized AM after 100 
iterations. 


III. Results 

The multiresolution technique has been evaluated using a 
real data scan of the NIST Phantom Test Article A HD 
acquired on a SureScan™ vlOOO Explosive Detection System. 
A two dimensional Level 3 Haar disrete wavelet transform is 
used to represent each z-slice of the three dimensional image 
domain. The wavelet tree, is initialized to consist of 
approximation coefficients only. At iteration number 64, the 
coefficients are back projected to the image space, voxel values 
across z-slices are summed up and the pixels whose values 
were larger than 0.1 times the maximum of the summed image 
were chosen to expand one level. Then, at iteration number 
128, the same procedure is applied with the same factor to 
expand one level further, and the last expansion is done at 
iteration number 256. Figure[^ shows objective function values 
versus time for unregularized alternating minimization algo¬ 
rithm (AM) Q and unregularized wavelet AM represented in 
this paper. AM algorithm has been run for 100 iterations while 
Wavelet AM has been run for 300 iterations. Figures and 

show image slices reconstructed from two algorithms at the 
same objective function value level. The difference between 
these two images (unregularized AM image subtracted from 
wavelet AM image) is shown in Figure It is important to 
note that even though two images are at the same objective 
function value level, the image reconstructed using wavelet 
AM has sharper edges. 
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Fig. 3: Image reconstructed with wavelet AM after 300 itera¬ 
tions. 

IV. Conclusion 

A fast, iterative, and adaptive algorithm for x-ray imaging 
was formulated and presented by using alternating minimiza¬ 
tion framework. The algorithm is guaranteed to decrease at 
each iteration and adaptive wavelet tree structure provides bet¬ 
ter utilization of computations. In other words, more computa¬ 
tions are used for the regions with high frequency components 
like edges while less are used for smoother areas. The wavelet 
tree expansion used to reconstruct the image shown in the 
results section is one of many possible methods to perform it. 
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Fig. 4: Difference image, unregularized AM image subtracted 
from wavelet AM image. 


Different ways to expand the tree will be investigated in the 
future. Different scale levels of discrete wavelet transform, 
different wavelet types and exploration of regularization are 
other parts to be explored later. Furthermore, this method can 
be combined with other acceleration methods like ordered sub¬ 
sets O- Preliminary studies combining ordered subsets and 
wavelet AM showed promising results and will be investigated 
further. 
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